/*
 * @Author: chyang0822 270917365@qq.com
 * @Date: 2024-10-14 14:06:21
 * @LastEditors: chyang0822 270917365@qq.com
 * @LastEditTime: 2024-10-18 17:37:09
 * @FilePath: /numpy_c_api/04-Py_SplineFilter1D/main.c
 * @Description: 
 * 
 * Copyright (c) 2024 by ${git_name_email}, All Rights Reserved. 
 */
// numpy_demo.cpp 
#include <Python.h> 
// #include <iostream>
#include <numpy/arrayobject.h>
// #include "ni_interpolation.h"


/* Converts a Python array-like object into a behaved input array. */
static int
NI_ObjectToInputArray(PyObject *object, PyArrayObject **array)
{
    int flags = NPY_ARRAY_ALIGNED | NPY_ARRAY_NOTSWAPPED;
    *array = (PyArrayObject *)PyArray_CheckFromAny(object, NULL, 0, 0, flags,
                                                   NULL);
    return *array != NULL;
}

/* Like NI_ObjectToInputArray, but with special handling for Py_None. */
static int
NI_ObjectToOptionalInputArray(PyObject *object, PyArrayObject **array)
{
    if (object == Py_None) {
        *array = NULL;
        return 1;
    }
    return NI_ObjectToInputArray(object, array);
}
/* Converts a Python array-like object into a behaved output array. */
static int
NI_ObjectToOutputArray(PyObject *object, PyArrayObject **array)
{
    int flags = NPY_ARRAY_BEHAVED_NS | NPY_ARRAY_WRITEBACKIFCOPY;
    /*
     * This would also be caught by the PyArray_CheckFromAny call, but
     * we check it explicitly here to provide a saner error message.
     */
    if (PyArray_Check(object) &&
            !PyArray_ISWRITEABLE((PyArrayObject *)object)) {
        PyErr_SetString(PyExc_ValueError, "output array is read-only.");
    return 0;
    }
    /*
     * If the input array is not aligned or is byteswapped, this call
     * will create a new aligned, native byte order array, and copy the
     * contents of object into it. For an output array, the copy is
     * unnecessary, so this could be optimized. It is very easy to not
     * do NPY_ARRAY_UPDATEIFCOPY right, so we let NumPy do it for us
     * and pay the performance price.
     */
    *array = (PyArrayObject *)PyArray_CheckFromAny(object, NULL, 0, 0, flags,
                                                   NULL);
    return *array != NULL;
}

/* The different boundary conditions. The mirror condition is not used
     by the python code, but C code is kept around in case we might wish
     to add it. */
// typedef enum {
//     NI_EXTEND_FIRST = 0,
//     NI_EXTEND_NEAREST = 0,
//     NI_EXTEND_WRAP = 1,
//     NI_EXTEND_REFLECT = 2,
//     NI_EXTEND_MIRROR = 3,
//     NI_EXTEND_CONSTANT = 4,
//     NI_EXTEND_GRID_WRAP = 5,
//     NI_EXTEND_GRID_CONSTANT = 6,
//     NI_EXTEND_LAST = NI_EXTEND_GRID_WRAP,
//     NI_EXTEND_DEFAULT = NI_EXTEND_MIRROR
// } NI_ExtendMode;

#include "ni_interpolation.h"
#include "ni_support.h"
#include "ni_splines.h"
#include <stdlib.h>
#include <math.h>

// int test(PyArrayObject *input, PyArrayObject* zoom_ar)
// {
//     npy_intp siz1 = PyArray_SIZE(input);
//     return 0;
// }


/* map a coordinate outside the borders, according to the requested
     boundary condition: */
static double
map_coordinate(double in, npy_intp len, int mode)
{
    if (in < 0) {
        switch (mode) {
        case NI_EXTEND_MIRROR:
            if (len <= 1) {
                in = 0;
            } else {
                npy_intp sz2 = 2 * len - 2;
                in = sz2 * (npy_intp)(-in / sz2) + in;
                in = in <= 1 - len ? in + sz2 : -in;
            }
            break;
        case NI_EXTEND_REFLECT:
            if (len <= 1) {
                in = 0;
            } else {
                npy_intp sz2 = 2 * len;
                if (in < -sz2)
                    in = sz2 * (npy_intp)(-in / sz2) + in;
                // -1e-15 check to avoid possibility that: (-in - 1) == -1
                in = in < -len ? in + sz2 : (in > -1e-15 ? 1e-15 : -in) - 1;
            }
            break;
        case NI_EXTEND_WRAP:
            if (len <= 1) {
                in = 0;
            } else {
                npy_intp sz = len - 1;
                // Integer division of -in/sz gives (-in mod sz)
                // Note that 'in' is negative
                in += sz * ((npy_intp)(-in / sz) + 1);
            }
            break;
        case NI_EXTEND_GRID_WRAP:
            if (len <= 1) {
                in = 0;
            } else {
                // in = len - 1 + fmod(in + 1, len);
                in += len * ((npy_intp)((-1 - in) / len) + 1);
            }
            break;
        case NI_EXTEND_NEAREST:
            in = 0;
            break;
        case NI_EXTEND_CONSTANT:
            in = -1;
            break;
        }
    } else if (in > len-1) {
        switch (mode) {
        case NI_EXTEND_MIRROR:
            if (len <= 1) {
                in = 0;
            } else {
                npy_intp sz2 = 2 * len - 2;
                in -= sz2 * (npy_intp)(in / sz2);
                if (in >= len)
                    in = sz2 - in;
            }
            break;
        case NI_EXTEND_REFLECT:
            if (len <= 1) {
                in = 0;
            } else {
                npy_intp sz2 = 2 * len;
                in -= sz2 * (npy_intp)(in / sz2);
                if (in >= len)
                    in = sz2 - in - 1;
            }
            break;
        case NI_EXTEND_WRAP:
            if (len <= 1) {
                in = 0;
            } else {
                npy_intp sz = len - 1;
                in -= sz * (npy_intp)(in / sz);
            }
            break;
        case NI_EXTEND_GRID_WRAP:
            if (len <= 1) {
                in = 0;
            } else {
                in -= len * (npy_intp)(in / len);
            }
            break;
        case NI_EXTEND_NEAREST:
            in = len - 1;
            break;
        case NI_EXTEND_CONSTANT:
            in = -1;
            break;
        }
    }

    return in;
}

#define BUFFER_SIZE 256000
#define TOLERANCE 1e-15



/* copy row of coordinate array from location at _p to _coor */
#define CASE_MAP_COORDINATES(_TYPE, _type, _p, _coor, _rank, _stride) \
case _TYPE:                                                           \
{                                                                     \
    npy_intp _hh;                                                     \
    for (_hh = 0; _hh < _rank; ++_hh) {                               \
        _coor[_hh] = *(_type *)_p;                                    \
        _p += _stride;                                                \
    }                                                                 \
}                                                                     \
break

#define CASE_INTERP_COEFF(_TYPE, _type, _coeff, _pi, _idx) \
case _TYPE:                                                \
    _coeff = *(_type *)(_pi + _idx);                       \
    break

#define CASE_INTERP_OUT(_TYPE, _type, _po, _t) \
case _TYPE:                                    \
    *(_type *)_po = (_type)_t;                 \
    break

#define CASE_INTERP_OUT_UINT(_TYPE, _type, _po, _t)  \
case NPY_##_TYPE:                                    \
    _t = _t > 0 ? _t + 0.5 : 0;                      \
    _t = _t > NPY_MAX_##_TYPE ? NPY_MAX_##_TYPE : t; \
    _t = _t < 0 ? 0 : t;                             \
    *(_type *)_po = (_type)_t;                       \
    break

#define CASE_INTERP_OUT_INT(_TYPE, _type, _po, _t)   \
case NPY_##_TYPE:                                    \
    _t = _t > 0 ? _t + 0.5 : _t - 0.5;               \
    _t = _t > NPY_MAX_##_TYPE ? NPY_MAX_##_TYPE : t; \
    _t = _t < NPY_MIN_##_TYPE ? NPY_MIN_##_TYPE : t; \
    *(_type *)_po = (_type)_t;                       \
    break



void display(PyArrayObject *mat, int dims)
{
    PyArrayObject *mat_3;
    PyArray_OutputConverter(mat, &mat_3);
    npy_intp *shape = PyArray_SHAPE(mat_3);
    double *array_3 = (double*)PyArray_DATA(mat_3);

    printf("numpy result:\n");
    if(dims == 2)
    {
        for (int i = 0; i < shape[0]; i++) {
        for (int j = 0; j < shape[1]; j++) {
            printf("%f\t", array_3[i*shape[1] + j]);
        }
        printf("\n");
        }
    }
    if (dims == 1)
    {
        for (int i = 0; i < shape[0]; i++) {
            printf("%f\t", array_3[i]);
        }
        printf("\n");
    }
    
}

int main(int argc, char *argv[])
{
    wchar_t *program = Py_DecodeLocale(argv[0], NULL);
    if (program == NULL) {
        fprintf(stderr, "Fatal error: cannot decode argv[0]\n");
        exit(1);
    }
    Py_SetProgramName(program);  /* optional but recommended */
    Py_Initialize();
    import_array();/* load numpy api */
    double array_input[3][3] = { { 2,5,5 },{ 5,6,8} ,{ 3,2,6}};
    npy_intp dims_input[] = { 3, 3 };
    PyArrayObject *mat_input = PyArray_SimpleNewFromData(2, dims_input, NPY_DOUBLE, array_input);

    double array_output[3][3] = {{ 0,0,0},{ 0,0,0} ,{ 0,0,0}};
    npy_intp dims_output[] = { 3, 3};
    PyArrayObject *mat_output = PyArray_SimpleNewFromData(2, dims_output, NPY_DOUBLE, array_output);

    double array_output2[3][3] = {{ 0,0,0},{ 0,0,0} ,{ 0,0,0}};
    npy_intp dims_output2[] = { 3, 3};
    PyArrayObject *mat_output2 = PyArray_SimpleNewFromData(2, dims_output2, NPY_DOUBLE, array_output2);

    display(mat_input, 2);
    
    NI_SplineFilter1D(mat_input, 3, 0, 0, mat_output);
    NI_SplineFilter1D(mat_output, 3, 1, 0, mat_output2);
    display(mat_output2, 2);

    // double array_input_3d[2][2][2] = { 2,5, 5,5, 6,8 , 3,2};
    // npy_intp dims_input_3d[] = { 2, 2, 2 };
    // PyArrayObject *mat_input_3d = PyArray_SimpleNewFromData(3, dims_input_3d, NPY_DOUBLE, array_input_3d);

    // double array_output_3d[4][4][4] = {0};
    // npy_intp dims_output_3d[] = { 6, 3};
    // PyArrayObject *mat_output_3d = PyArray_SimpleNewFromData(3, dims_output_3d, NPY_DOUBLE, array_output_3d);

    
    // double array_zoom_3d[] = {0.5, 0.5, 0.5};
    // npy_intp dims_zoom_3d = 3;
    // PyArrayObject *mat_zoom_3d = PyArray_SimpleNewFromData(1, &dims_zoom_3d, NPY_DOUBLE, array_zoom_3d);
    // display(mat_zoom, 1);

    
    // npy_intp siz1 = PyArray_SIZE(mat_output);
    // display(mat_input, 2);
    
    // NI_ZoomShift(mat_input_3d, mat_zoom_3d, Py_None, mat_output_3d, 3, (NI_ExtendMode)0, 0.0, 0, 1);
    // display(mat_output, 2);

    
    

    return 0;
}